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Abstract. In this paper we propose a blind deconvolution method which applies to 
data perturbed by Poisson noise. The objective function is a generalized Kullback-Leibler 
(KL) divergence, depending on both the unknown object and unknown point spread 
function (PSF), without the addition of regularization terms; constrained minimization, 
with suitable convex constraints on both unknowns, is considered. The problem is non- 
convex and we propose to solve it by means of an inexact alternating minimization 
method, whose global convergence to stationary points of the objective function has 
been recently proved in a general setting. The method is iterative and each iteration, 
also called outer iteration, consists of alternating an update of the object and the PSF by 
means of fixed numbers of iterations, also called inner iterations, of the scaled gradient 
projection (SGP) method. Therefore the method is similar to other proposed methods 
based on the Richardson-Lucy (RL) algorithm, with SGP replacing RL. The use of SGP 
has two advantages: first, it allows to prove global convergence of the blind method; 
secondly, it allows the introduction of different constraints on the object and the PSF. 
The specific constraint on the PSF, besides non-negativity and normalization, is an 
upper bound derived from the so-called Strehl ratio (SR), which is the ratio between the 
peak value of an aberrated versus a perfect wavefront. Therefore a typical application, 
but not the unique one, is to the imaging of modern telescopes equipped with adaptive 
optics systems for partial correction of the aberrations due to atmospheric turbulence. 
In the paper we describe in detail the algorithm and we recall the results leading to its 
convergence. Moreover we illustrate its effectiveness by means of numerical experiments 
whose results indicate that the method, pushed to convergence, is very promising in the 
reconstruction of non-dense stellar clusters. The case of more complex astronomical 
targets is also considered, but in this case regularization by early stopping of the 
outer iterations is required. However the proposed method, based on SGP, allows 
the generalization to the case of differentiablc regularization terms added to the KL 
divergence, even if this generalization is outside the scope of this paper. 
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1. Introduction 

Blind deconvolution is the problem of image deblurring when the blur is unknown and, 
in general, is investigated by assuming a space-invariant model; in such a case the naive 
problem formulation is to solve the equation 

g = h*f , 

where g is the detected image and (/, h) are respectively the unknown object and the 
unknown point spread function (PSF), while * denotes convolution. It is obvious that 
this problem is extremely undetermined and that there is an infinite set of pairs solving 
the equation. Among them also the trivial solution (f = g , h = S), where <5 denotes the 
usual delta function. Therefore the problem must be reformulated by introducing as far 
as possible all available a priori information on both the object and the PSF. 

Blind deconvolution is the subject of a wide literature and we do not try to give 
a thorough account for that. Indeed the different approaches concern specific classes 
of images and PSFs. For instance approaches applicable to natural images may not be 
suitable in microscopy; approaches developed for motion blur are not applicable to other 
classes of blur, and so on. As concerns natural images we only mention a recent paper 
[32] which contains a critical analysis as well as several relevant references. 

In this paper we focus on astronomical imaging by assuming that an adaptive 
optics (AO) system is used to compensate for atmospheric blur and that a parameter 
characteristic of this correction, the so-called Strehl ratio (SR), is approximately known. 
We recall that SR is the ratio of peak diffraction intensity of an aberrated versus perfect 
waveform. In the case of AO images this parameter can be estimated by the astronomers 
during the observation and provided with an error of few percent (about 4-5 %). Since 
this information provides an upper bound on the maximum value of the PSF, it can be 
used to exclude the trivial solution mentioned above and corresponding to the pair (g, S). 

The approach we propose applies to astronomical imaging if noise is dominated by 
photon counting and therefore the data are realizations of Poisson processes, even if 
approaches based on regularized least-square methods are also available (see for instance 
[29]). In the Poisson case several iterative methods have been already investigated, which 
consist of alternating updates of the object and PSF by means of Richardson- Lucy (RL) 
iterations [2B1 1231 EDI S21 1201 121] , or accelerated RL iterations ^Q\. 

In [28] one iteration of the algorithm consists of updating both the object and the PSF 
by means of one RL iteration. This algorithm was investigated, in a different context, by 
Lee and Seung [31] but their convergence proof is incomplete, since only the monotonic 
decrease of the objective function is shown while, for a general descent method to be 
convergent, strongest Armijo-like decreasing conditions have to be verified [35]. The 
other approaches could be classified as methods of inexact alternating minimization since 
they use alternately a number of RL iterations on both the object and the PSF (remark, 
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however, that the optimization problem underlying these approaches is not explicitly 
mentioned by the authors). Their convergence is not proved if RL, or the acceleration of 
RL proposed in [10] . are the algorithms used for inexact optimization. 

In a recent paper [2] , in the context of non- negative matrix factorization, convergence 
of inexact alternating minimization is proved if the iterative algorithm used for the inner 
iterations satisfies suitable conditions, which are satisfied by the scaled gradient projection 
method (SGP) [11] proposed for constrained minimization of convex differentiable 
functions. Therefore in this paper we utilize these results to propose a convergent 
blind deconvolution approach applicable to the reconstruction of astronomical adaptive- 
optics (AO) corrected images. We remark that the advantage of SGP is not only a 
fast convergence, if a suitable scaling of the gradient is used, but also the possibility of 
introducing suitable convex constraints on the solution. The practical limitation is that 
the projection operator onto the convex set defined by the constraints should be easily 
computable. This is the case of box and equality constraints and the constraints we 
introduce on the PSF just belong to this class. Remark that, in the case of Gaussian 
noise, a similar situation is achievable if the projected Landweber method is used (for 
an application to seismology see [5]) or, more precisely, the accelerated version provided 
by the application of gradient projection methods [3]. Indeed, in this case the conditions 
required in [T3] for convergence are satisfied. 

The structure of our approach is similar to that of the previously mentioned methods 
based on the RL algorithm, the difference being of course that RL is replaced by SGP 
with different constraints on / and h: in the case of the object we only consider non- 
negativity while in the case of the PSF we consider both non-negativity and an upper 
bound provided by the knowledge of the SR, as well as the normalization condition which 
must be satisfied by the PSF. We point out that the relevance of the use of the SR 
constraint for blind deconvolution was first pointed out by Desidera & Carbillet [21] and 
this paper intends to use it in a proper mathematical context. 

Thanks to [H] the convergence is assured if we use a fixed number of SGP iterations 
for updating the object and the PSF; the number of iterations may be different in the two 
cases (for the denomination asymmetric iterative blind deconvolution, see [Hi]). Since the 
problem is non-convex, the limit of the iteration may depend on the choice of the initial 
step and possibly on the numbers of internal iterations. The convergence result does 
not assure that the limit is a sensible solution of the problem, since we do not introduce 
regularization in our approach. A comment on this point is required. 

In the case of deconvolution it is well known that the minimizers of the discrepancy 
function for Poisson data, the generalized Kullback-Leibler (KL) divergence, are sparse 
objects, i.e. they consist of bright spots over a black background; it is the so-called night- 
sky [1] or checker-board [31] effect. As a result these minimizers can be sensible solutions 
in the case, for instance, of the deconvolution of images of not dense star clusters by a 
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given PSF and this result is confirmed by a wide set of numerical experiments. On the 
other hand, if the data is the image of a star cluster and we deconvolve it using a sparse 
object with points correctly located at the positions of the stars, we may expect that the 
result is a satisfactory reconstruction of the PSF; this reconstruction should improve as 
the reconstructed image used in the deconvolution improves. Therefore, using a suitable 
strategy in the choice of the number of inner iterations, we can expect sensible results in 
the case of stellar objects, by pushing to convergence the outer iterations. This argument 
is supported by our numerical experiments. 

The situation is different in the case of diffuse or complex objects. In this case we 
can believe that the semi-convergent behaviour of RL or SGP or similar methods implies 
a similar behaviour for the outer iterations of the blind method; and this is just what we 
find in our tests. An alternative is obviously the introduction of regularization by adding 
suitable penalty terms to the KL divergence. However there are two main problems: 
the first is the selection of a suitable regularization for the astronomical object to be 
reconstructed and of a suitable regularization for the PSF to be reconstructed; the second 
is the selection of a suitable rule for estimating the two regularization parameters. We 
do not know generally accepted solutions for both problems and therefore early stopping 
of the iterations is still the easiest approach to regularization. We only remark that the 
proposed method can be easily generalized to the case of differentiable penalties thanks to 
the generalization of SGP to the regularized proposed in several papers [4"3 | PfO l fT3] . 

The paper is organized as follows. In Sect. 2 we formulate the blind deconvolution 
problem as a constrained minimization of the generalized KL divergence as follows from 
a maximum-likelihood approach to Poisson data deblurring [6]. In Sect. 3 we summarize 
the results on the inexact alternating minimization problem proved in [H] as well the 
main features of the SGP method proposed in [TT]. In both cases we also provide the 
algorithms used in this paper. Finally in Sect. 4 we describe our numerical experiments 
with a particular focus on the case of astronomical objects consisting of small star 
clusters, represented in our simulations by point-wise objects. In these cases we observe a 
remarkable convergence of the reconstructed PSF to that used in image simulations. We 
also attempt an accurate analysis of the artifacts generated in the reconstructed images 
since an understanding of their structure may be important in the practical applications 
of our method. Sect. 5 is devoted to a few conclusions and possible extensions. 

2. Problem setting 

Following [3E], we assume that the observed image g can be modeled as the sum of two 
terms g = g pe + r . The first, g , is the number of photo-electrons due to object and 
background emission and is a realization of a Poisson random variable with expected value 
g = h * f + b, where / is the original object, h is the PSF of the acquisition system and 
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b is the background term, while r represents the read-out noise (RON). Here and in the 
following we denote by bold letters N x N arrays whose pixels are indexed by a multi- 
index i = (21,22), i 6 5. For simplicity we assume that the background is constant and 
known. As concerns the RON, it is a realization of a Gaussian additive random variable 
with a known variance a 2 . According to Snyder et al. [39], it can be approximated by a 
Poisson process with mean and variance being the same as a 2 if the constant term a 2 is 
added to g. If we add a 2 also to the background and if, with an abuse of notation, we 
denote again as g and b the modified image and background, then we can conclude that 

g ~ Poiss(/i * f + b) . 

As concerns the PSF, we assume that it is normalized to unit volume 

ies 

and that its maximum value, denoted by s, is known 

max hi = s. (2) 

ies 

The upper bound s can be obtained by computing the diffraction-limited PSF of the 
considered telescope and multiplying its peak value by the SR value provided by the 
astronomers, as discussed in the Introduction. 

The blind deconvolution problem consists in finding an approximation of both / and 
h, given g, b and s. To this purpose we consider a maximum-likelihood approach to the 
problem of image deconvolution. Since the maximization of the likelihood, which depends 
on the unknown object and PSF, is equivalent to the minimization of a generalized KL 
divergence, we propose to estimate these approximations by minimizing this function (see 
the comments in the Introduction concerning regularization) while taking into account all 
the available information, i.e. the non-negativity of both the PSF and the original object 
and the constraints (JTJ)- (J2J) . The resulting optimization problem is the following 

min KL(g,h* f + b) (3) 
s.t. / > ; 0<h<s ,^/ii = l 

where KL denotes the KL divergence of h * f + b from g 

KL(g,h*f + b) = J2\si l °S (h fl + b . + ( h *f)i + b i-9i\ ■ (4) 

ieS I \ * J)i + i J 

Problem ([3]) is convex if restricted to / or h only, but is in general non-convex with 
respect to the pair (/, h), thus leading to the possible presence of several local minima. 
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Indeed, the gradient and Hessian of the objective function in are given by 



V f KL(g, h*f + b) = l-H 
V h KL(g,h*f + b) = F T l 
V 2 KL(g, h*f + b) = 



T 



h*f+b 



71T 



h*f+b 



( 



H diag 



(h*f+by 



H 



iJ T diag 



\ 



F T ^g (jf^rw) H - K(h, ff i^diag 



F-K(hJ) > 
F 

I 



-b) 



where H and F are the Block Circulant with Circulant Blocks (BCCB) matrices associated 
to the convolution, i.e. h* f = Hf = Fh, while K(h, /) is the Block Hankel with Hankel 
Blocks (BHHB) matrix whose last row is the vector h ^ +b (see [27J Chapter 4] for a survey 
on structured matrices). Here the ratios and the squares are computed element-wise, and 
1 is a column vector with all entries equal to 1. Even if the diagonal blocks of the Hessian 
are symmetric positive semi-definite, the whole matrix is difficult to analyze and compute. 



3. Alternating Minimization 

Despite the complexity of the Hessian, the constraints have a simple, separable structure, 
which can be exploited by adopting an Alternating Minimization (AM) algorithm for the 
solution of the non-convex problem (JHJ)- . More precisely, the AM algorithms can be 
applied to any problem of the form 

min J(x) (5) 
s.t. x e fii x f2 2 x ... x Vt m C W 1 

where, for allz = 1, m, f2j is a closed and convex subset of R ni with n\ + ... + n m = n 
and any vector in the feasible set can be partitioned into vector components as x = 
(xi, x 2 , x m ) Xi G fij. Clearly, the blind deconvolution problem ()3]) is a special case 
of ([5]) with rn — 2, X\ — f and x<i = h. 

The basic idea of AM is the cyclic minimization of the objective function with 
respect to one variable, updating its value for the next optimization steps: in particular, 
AM is often referred to as the Nonlinear Gauss-Seidel ( GS) method, where the iterate 
x {k+i) = (^xf +l \ ... : Xm +l ^) is computed such that for i = l,...,m the block of variables 
cc 4 is a solution of the sub-problem 

min J{x^ 1 \...,xlf\y,x%...,xM). (6) 
s.t. y e fit 

This kind of approach has been widely studied in the literature [U [91 |25j |26j |33j |JT] and 
we recall two important facts about it: 
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• for m — 2 it has been proved in [2TH Corollary 2] that the limit points of the sequence 
{x^} defined in ([6]) are stationary for problem ([5]) even in the non-convex case; 

• for m > 3 the convergence of the nonlinear GS method (jH]) to a solution of §5§ is 
not guaranteed, without additional convexity assumptions on the objective function 
J: indeed, in [36], Powell devises a counterexample with m = 3 where all the limit 
points of the sequence generated by the nonlinear GS method are not stationary for 
the problem ([5]). Some convergence results are proved for example in P [9J EHJ [33] 
under suitable strict convexity assumptions. 

All the convergence results mentioned above, even in the case m = 2, are proved when the 
iterates are updated by an exact solution of the partial minimization problem ([6]), which is 
often impractical or too costly to compute. Indeed, many practical AM algorithms, which 
are also referred to as Block Coordinate Descent methods, are obtained by applying an 
iterative minimization method to approximately solve fl6]). In this case, the convergence 
properties of the alternating scheme also depend on the features of the inner solver. A 
detailed analysis of the Block Coordinate Descent algorithms in the unconstrained case 
is proposed in [25] , where the authors devise some convergence conditions not necessarily 
related to the convexity of the objective function. 

In this paper, we follow the approach in [14], where the partial minimization over 
each variable (Q is performed inexactly by means of a fixed number of Scaled Gradient 
Projection (SGP) steps [TT] . 

Choose a feasible starting point x^ and a positive integer L > 1 
For fc = 0,1,2, ... 

(7) 

For % = 1 , . . . , m v ' 

Compute x[ k+1 ^ by applying n\ k < L SGP iterations to Q) 

A representation of the scheme ([7]) applied to the blind deconvolution problem is 
given in Algorithm [1] each main cycle consists of two successive deconvolution steps to 
update the current estimates of the object and PSF tv k \ For sake of completeness, 
we report the convergence result shown in Theorem 4.2 of [H] for our particular case. 

Theorem 3.1 Every limit point of the sequence {f^ k \h^) generated by Algorithm^ is 
a stationary point for problem (G|). 

As far as we know, convergence results stronger than that given in this Theorem for 
first-order methods applied to a general non-convex problem do not exist in the literature. 
The main difficulty in this kind of problems is to prove the existence of convergent sub- 
sequences. However, in our specific case, thanks to the Strehl constraint, the sequence 
of the PSFs generated by our approach is bounded. Moreover, as concerns object 
reconstruction one could introduce a constraint on the flux (£i norm) of the reconstructed 
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objects since this parameter can be derived from the data. However, as follows from 
our numerical experience this constraint is practically assured by SGP (we recall that 
it is exactly assured by RL in the case of zero background) and therefore we do not 
introduce it in order to reduce the computational cost. It follows that also the sequence 
of the reconstructed objects is bounded. We conclude that the existence of convergent 
sub-sequences is assured. We can add that we remarked a convergent behaviour of all the 
sequences obtained in our numerical experiments. 

Finally, we point out that the convergence result holds for any number of inner 
SGP iterations. The key point of such theoretical analysis is the sufficient decrease of 
the objective function, which is enforced at each SGP iteration by means of an Armijo 
backtracking loop. Since the objective function (j4]) is not convex with respect to the 
couple (/, h), the presence of multiple potential stationary points makes any limit point 
dependent on both the initial guess and the chosen inner iteration numbers on the image 
(nf) and the PSF (n h ). 



Algorithm 1 Cyclic Scaled Gradient Projection (CSGP) method 

Choose the starting point f^°\ and the inner iterations numbers nf,nh > 1. 

For k = 0, 1, 2, ... do the following steps: 

Step 1. Compute with n/ SGP iterations applied to 

min KL(g,h {k) * f + b) (8) 

s.t. / > 

starting from the point 
Step 2. Compute h^ k+1 ^ with nh SGP iterations applied to 

min KL{g,h* f k+1) + b) (9) 

s.t. 0<h<s ,Y^i = l 

ies 

starting from the point h^ k \ 

End 



We stress the fact that the main strength of Algorithm [T] for blind deconvolution with 
respect to the more standard AM approach described for example in |16], is that it allows 
an inexact solution of the inner sub-problems ©-Q while preserving the theoretical 
convergence properties. Since the proposed method is essentially based on SGP, we recall 
its main features in the following sub-section. 
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3.1. The Scaled Gradient Projection method 

The SGP algorithm is a first-order method which applies to any optimization problem of 
the form 

min J(x) (10) 

where J{x) is a continuously differentiable function and Q is a convex set. Each SGP 
iteration is based on the feasible descent direction defined as 

d {k) = Pn,D-^ (k) - a k D k VJ(x^)) - 

where a k is a scalar step-size parameter, D k is a diagonal matrix with positive diagonal 
entries and P n D -i(-) is the projection onto f2 associated to the norm induced by D^ 1 , i.e. 

P q,dA x ) = argmiri(a; - y) T D^\x - y). (11) 

The new point is computed along the direction d^ k ' as follows 
x (k+i) = x (k) + Xkd (k) 

where X k is a step-length parameter to be chosen such that the monotone Armijo condition 

J(x {k) + \ k d {k) ) < J(x {k) ) + (3\ k VJ{x^) T d {k) (12) 

is satisfied for a fixed value of the parameter G (0, 1), in order to guarantee the sufficient 
decrease of the objective function. In practice, X k is computed by a standard backtracking 
condition as \ k = m , where 9 G (0, 1) and m is the smallest integer such that ( fl"2l) is 
satisfied. 

The convergence of the SGP scheme, which is outlined in Algorithm [2], can be proved 
when the step-size a k and the diagonal entries of D k are bounded above and away from 
zero, i.e. a k e [a min , a max ] with < a min < a max and D k is chosen in the set V 
of diagonal matrices whose diagonal entries have values between Li and L 2 , for given 
thresholds < L\ < L 2 . 

The SGP algorithm has been recently applied in several image restoration problems 
(see e.g. [31 H21 HH])- Under standard assumptions, it can be proved [11] that the SGP 
algorithm is well defined and any limit point of the sequence {x^} is a stationary point 
of (fTUj) ; if, in addition, J(x) is convex, any limit point is a minimum point. It is worth 
stressing that the convergence result holds for any choice of the scaling matrix D k G T> 
and the step-length a k G [a m i n ,a max ]: this freedom of choice can be exploited in order 
to improve the convergence speed. Indeed, it is well known that gradient methods can 
be significantly accelerated by a clever choice of the step-length parameter a k . one of the 
more effective strategies are the Barzilai-Borwein (BB) rules proposed firstly in [2] for 
quadratic unconstrained programming and then developed and analyzed for more general 
problems (see [T71[19l[2HI33] and reference therein). The BB rules can be considered a very 



Blind deconvolution in post- adaptive- optics astronomical imaging 



10 



Algorithm 2 Scaled gradient projection (SGP) method 

Choose the starting point x^ > and set the parameters j3, 9 G (0, 1), < a min < a max . 

For k — 0, 1, 2, ... do the following steps: 

Step 1. Choose the parameter a k G [a m m, a m «i] and the scaling matrix D k G V; 
Step 2. Compute the descent direction: 

rfW = P QD -i(xW - a k D k VJ(xW)) - x^; 
Step 3. Backtracking loop: compute the smallest integer m such that (fT2l is 

satisfied with \ k = 9 m ; 
Step 4. Set x^ = x^ + \ k d {k) . 

End 



cheap way to capture the second order information enforcing a quasi-Newton property. 
In our algorithm we adopt the scaled versions of the BB rules proposed in [11] , which are 
given by 

where s^" -1 ) = a?( fe ) — and = VJ(a?( fc )) — V J(x^ k ~^). Based on the previous 

formulas, we define the values a k ,a k G [a m i n ,a max ] in the following way 

IF s^-^D^z^ < THEN 
a k 1] = min {10 • a fc _i, a max }; 



ELSE 



a k ' = mm {a max , max jaw, a fc 

ENDIF 

IF D k Z^ k ^ < THEN 

„(2) 



-mi- 

T - 



a\ = min {10 • a fc _i, awj; 

ELSE 

a^. 2) = min {a max , max|a min ,, ajjf^}}; 

ENDIF 

From the numerical experience, the best performances are obtained by an alternation of 
the two BB formulas: thus, following [21], we choose the following criterion for computing 

a k 

IF k < 20 THEN 



(1) 

a k = min a) ; (14) 

j=max{l,k+l-M a },...,k J 



ELSE IF a k 2 '/a k < T k THEN 

Set a k as in (fT4l 
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r k+ i = 0.9 ■ r k ; 



ELSE 



a k = a k 



(i). 



ENDIF 



where M a is a prefixed positive integer and n G (0, 1). Our choice of taking the value 
defined in ( fT4|) for the first 20 iterations leads to a more stable behaviour and, in some 
cases, also to a slight improvement of the reconstruction accuracy [37] . 

The other crucial ingredient for the practical performances of SGP is the choice of 
the scaling matrix: in our case, taking into account the objective functions of ([S]) and fl9]), 
we adopt the scaling suggested by the RL algorithm 



as suggested also in [IT]. Obviously, in the inner iterations of flS}, a; is replaced by / 
while, in the inner iterations of ([9]), a? is replaced by h. As concerns the choice of the 
bounds (Li,L 2 ), at the beginning of each inner subproblem we perform one step of the 
RL method and tune the parameters according to the min/max positive values ?/mm/?/max 
of the resulting image according to the rule 



ENDIF 

3.2. Computing the projections 

Since the minimization steps (|8]) and (Q involve different constraints, corresponding to 
two convex sets Qi and Q 2 , respectively, we have to account for two different algorithms to 
compute the projections P n D -i and P n D -i. In the alternating procedure of Algorithm 
[H when SGP is applied to problem ([8]), the projection consists of a simple component 
thresholding, obtained by setting all the negative elements of the vector to be projected 
equal to zero. For the updating of the PSF, instead, we have to project on the constraints 
set of the problem ([9]), consisting of a single linear equality constraint, in addition to 
simple bounds (box constraints) on the variables. The resulting constrained optimization 
problem to be addressed is therefore 



D k = diag (min(L 2 , max(Li, ar fc ))) 



IF Z/max/ymin < 50 THEN 

Li = y m in/W; 
L 2 = 2/ma X ■ 10; 



ELSE 



L\ — ymin'i 
^2 = Umax] 



min -z D k z — z y 



(15) 



s.t. < z < s 
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where y = D~^ l (x^ — akD^V J(x^)). By introducing the Lagrangian penalty function, 
one can see that the orthogonal projection (fl5l) can be re-conducted to a root-finding 
problem of the piecewise linear monotonically non-decreasing function 

*(o = 5>^)- 1 = ' 

where £ is the Lagrangian multiplier of the equality constraint, 

*i(0 = mid(0 J (D fc ) ii ( ! / i + 0,5) 

and mid(ai, 02, (J3) is the component-wise operation that supplies the median of its three 
arguments. For solving this kind of problem we apply the secant-based method proposed 
in [18] (its Matlab implementation is given in the corresponding technical report download- 
able at the webpage http://www.maths.dundee.ac.uk/nasc/na-reports/NA216_RF.pdf), 
which is able to compute the projection very quickly and whose computational cost grows 
linearly in time with respect to the image size [ITJ §3.1]. 



4. Numerical experiments 

In this section we investigate the effectiveness of the proposed blind method by means of 
several numerical experiments. Since the blind problem formulated in is non-convex, 
several local minima may exist. Moreover, we know that any limit point (or the limit) 
of the proposed iteration is a stationary point of the problem. The limit depends, in 
general, on the numbers of inner iterations but also on the initialization of the outer 
iteration; therefore it is important to initialize the procedure with a sensible initial guess 
and we first discuss this point. 

As concerns the object we can use the standard initialization of the RL algorithm, 
namely a constant object with a flux coinciding with the flux of the image after background 
subtraction. The choice of the initial PSF is more important because, in the first step of 
the procedure, the image is deconvolved with this PSF. 

To this purpose we point out an important property of the PSF of a telescope: it is a 
band-limited function and, if the telescope consists of a circular mirror, the band, i.e. the 
support of its Fourier transform, is a disc with a radius proportional to the ratio between 
the diameter D of the telescope and the observation wavelength A. It is not easy to insert 
this property as a constraint on the PSF because the projection on the resulting set of 
constraints (including SR and normalization) is not easily computable. For this reason we 
do not consider this constraint in this paper. However we can try to force the estimated 
PSF to have this property using an initial PSF which is band-limited and satisfies the 
other constraints. 

The ideal PSF of the telescope is not suitable as initial guess because it does not 
satisfy the SR constraint. However one can consider, as suggested for instance in [TO] . 
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the autocorrelation of the ideal PSF, which has the same band. In our simulations, 
which assume in general a telescope of the 8m class and observations in H-band (see the 
beginning of the next section), the resolution, inversely proportional to the bandwidth, is 
about 50 mas (milliarcseconds) . Using an oversampling such that the pixel size is 15 mas, 
the diameter of the band in Fourier space is about 186 pixels (remark that it depends also 
on the number of pixels in the image). With these values, the autocorrelation is a good 
choice if SR > 0.46 (the value of s depends on both SR and the ratio D/X); for lower 
values of SR one can take the autocorrelation of the autocorrelation and so on, until the 
SR constraint is satisfied. This is the choice considered in our numerical experiments and, 
quite surprisingly, it seems that the algorithm, in spite of its high nonlinearity, preserves 
the band-limiting property satisfied by the initial guess. 

All the numerical experiments have been performed with a set of routines 
implemented by ourselves in Interactive Data Language (IDL). The codes of the algorithms 
presented and discussed in this paper are available under request. 

4-1. Image generation 

As mentioned in the Introduction the use of SR as a constraint on the PSF is first proposed 
in [21]. Therefore some of our numerical experiments coincide with some of the tests 
performed in that paper. In particular we use three of the AO-corrected PSFs (with 
SR equal to 0.67, 0.40 and 0.17, respectively), used by these authors and obtained by 
means of the Software Package CAOS [15] ; the parameters corresponding to these PSFs 
are given in [21]. We only specify that they correspond to a telescope with an effective 
diameter of 8.22 m and an observation wavelength of A = 1.65/im (H-band). For each 
PSF, the images are generated by assuming, as in [21] a time exposure of 1200 s, with a 
total transmission of 0.3. Moreover, a background of 13.5 mag arcsec -2 , corresponding to 
observations in H-band, is added to the blurred images (for the convenience of the reader 
we remark that it corresponds to 3.41 x 10 4 counts per pixel). The results are perturbed 
with Poisson noise and additive Gaussian noise with a = 10 e~/px. According to the 
approach proposed in [39] and discussed in Sect. 2, RON compensation is obtained in 
the deconvolution algorithms by adding the constant cr 2 = 100 to the images and the 
background. 

In our first experiment we also consider an example which is not related to AO 
imaging but is a simulation of HST image before COSTAR correction, since this image 
is frequently used in the testing of deconvolution methods. Obviously in such a case the 
ideal PSF must be computed, by taking into account the diameter of the Hubble telescope, 
about 2.4 m, and the assumed observation wavelength of about 0.55 /xm, and compared to 
the aberrated PSF in order to estimate the corresponding SR. Objects, PSFs and blurred 
images used in all our experiments are sized 256 x 256 pixels. 
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4-2. Point-wise objects 

We first report results on the following examples: 

• the binary system considered by Desidera & Carbillet [2T], in which the two 
components have the same magnitude 12 in H-band (corresponding approximately 
to 6.03 x 10 8 counts) with an angular separation of 285 mas (19 pixels), i.e. ~7 times 
larger than the diffraction limit (~40 mas); 

• a model of an open star cluster based on an image of the Pleiades, consisting of 9 
stars with magnitudes ranging from about 13 (i.e. about 2.32 x 10 8 counts) to 16 
(about 1.79 x 10 7 counts) in H-band and described in [7j; 

• a simulation of a star cluster, consisting of 470 light sources, as observed by the 
Hubble Space Telescope (HST) before COSTAR correction. For this case only, we 
do not use an AO-corrected PSF but the aberrated HST PSF, which corresponds to 
SR=0.09. These data can be obtained via anonymous ftp from 

|ftp: //ftp.stsci.edu/software/stsdas/testdata/restore/sims/star_cluster/ 

In Fig. 1 we show the images of the binary and of the star cluster in the case of a PSF 
with SR=0.67 as well as the HST image of a simulated star cluster. For all these examples 
we use 50 inner iterations on the object and 1 inner iteration for the PSF. This choice can 
be justified by the features of our blind problem discussed in the Introduction, since we 
need a sufficiently large number of SGP inner iterations for obtaining a nearly point-wise 
object. Moreover, with a few experiments on the binary, we verify that this choice is a 
good compromise which provides a sufficiently fast convergence for all cases. In a first 
instance we perform 300 outer iterations. 

As concerns the measure of the quality of the reconstructions, for the PSFs we use the 
relative r.m.s. error between the reconstructed PSF h and that used for image generation 
h, i.e. 

RMSE= l|/l T^ 12 , (16) 

where ||.||2 denotes the £2 norm. The same parameter is used for measuring the quality 
of the reconstruction of the HST star cluster. In the case of the binary and of the open 
star cluster we use a magnitude average relative error (MARE) defined as follows 

MARE= 1 -Y lm - Th * 1 , (17) 

where q is the number of stars and rrii, fhi are respectively the reconstructed and the true 
magnitudes. 

The results are shown in Table [T] and are consistent with the results reported in [21] 
but obtained with a sound mathematical approach, allowing investigation of the limit 
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Figure 1. Images of the binary and of the star cluster (upper panels) in the case of a 
PSF with SR—0.67. In the lower panel the image of the HST star cluster. 



for large number of iterations and generalization to regularized problems. The values 
of MARE estimated with our blind approach are certainly higher than those achievable 
if one deconvolves the data with the exact PSF (inverse crime) and given in the third 
column, but they are still quite small. Moreover the reconstruction of the PSFs is very 
satisfactory: the RMSEs of the initial PSFs are of the order of 30-50 %, while those of 
the reconstructed PSFs are of the order of few percents. A comparison between the true, 
initial and reconstructed PSFs is shown in Fig. 2. We must add that the reconstruction 
error is still decreasing after 300 iterations and therefore the minimum of the objective 
function is not yet reached. 

For investigating the limit of the algorithm we consider three other examples of 
binaries: a binary with magnitudes 12-12 and angular distance of 10 pixels and two 
binaries with magnitudes 12-16 and angular distances of 19 and 10 pixels respectively. 
For these four examples of binaries we generate images using the PSF with the highest 
SR, namely 0.67, and we compute 8000 outer iterations of the blind algorithm, using 
the fixed pair (nf,nh) = (50,1). In Table [2] we report the results obtained after 300, 
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Figure 2. First column: the PSFs used for image generation; second column: the 
PSFs used for initializing the blind algorithm (see the text for their computation); third 
column: the PSFs reconstructed by the blind algorithm. First row: AO-corrected PSF 
with SR=0.67; second row: AO-corrected PSF with SR=0.40; third row: AO-corrected 
PSF with SR=0.17; fourth row: HST PSF before COSTAR correction. 
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Table 1. Reconstruction errors for point-wise objects. In the first column we specify the 
object and in the second the value of the SR used for image generation; in the third and 
fourth the values of MARE (RMSE in the case of HST cluster) when SGP is used for 
image deconvolution with the exact and initial PSF respectively. In the fifth column the 
values of MARE (RMSE in the case of HST) obtained with 300 outer iterations by our 
blind approach, using the fixed pair (n/,n^) = (50, 1). Finally, in the last two columns 
the errors between the true and the initial PSF, followed by the errors between the true 
PSF and that provided by the blind approach. 



Image 


SR 


MARE 


MARE! 


MARE 2 


RMSEi 


RMSE 2 


Binary 


0.67 
0.40 
0.17 


1.86 x 10~ 5 
1.84 x 10~ 5 
2.36 x 10" 6 


1.12 x 10~ 2 
2.26 x 10" 2 
1.67 x 10~ 2 


1.44 x 10~ 3 
2.15 x 10~ 3 
1.99 x 10~ 3 


32 % 

54 % 

55 % 


1.8 % 

2.9 % 
3.3 % 


Cluster 


0.67 
0.40 
0.17 


2.10 x 10" 5 
4.43 x 10~ 5 
5.42 x 10~ 5 


1.07 x 10~ 2 
1.14 x lO" 2 
1.30 x lO" 1 


3.09 x 10" 4 
3.63 x 10" 4 
2.87 x 10~ 3 


32 % 

54 % 

55 % 


1.0 % 

1.1 % 

4.2 % 


HST 


0.09 


5.1 % 


25 % 


7.6 % 


47 % 


6.7 % 



4000 and 8000 outer iterations. We do not find a uniform behaviour: in some cases the 
errors are decreasing for increasing number of iterations while in others they are slightly 
increasing or do not have a monotone behaviour. In all cases the variations are small and 
the reconstruction errors on the PSF after 8000 iterations are quite small. However, in 
view of obtaining a sufficient accuracy, 300 outer iterations can be sufficient in all cases. 

Similar results are obtained in the case of the open star cluster. Moreover, since 
the previous examples are derived from the examples considered in [2T] and generated 
by assuming a very long observation time (hence a low noise level) we also generated 
an image of the binary 12-12, angular distance 19 pixels, assuming an integration time 
of 12 seconds, with a reduction by a factor 100 of the average number of photons. By 
performing 8000 iterations we still find convergence of the algorithm but the error on the 
PSF is now of the order of 1 % and can be reached after 300 iterations. On the other 
hand the value of MARE is quite satisfactory since it is of the order of 8 x 10 -4 . 

It should be interesting to find a way for establishing if the minima we find are the 
global ones or not, but, as it is known, global minimization is a very difficult problem. 
As a test, even if it does not provide a proof that the minima are the global ones, we 
compare the minimum values of the objective function, i.e. the KL divergence, with its 
values corresponding to the ground truths, i.e. the values obtained by substituting in Eq. 
01]) the objects and PSFs used for image generation. We find that the minimum values are 
of the order of 1.0 x 10 4 while the values corresponding to the ground truth are greater by 
about a factor of 3. We can only say that, if these values were smaller than our minimum 
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Table 2. Reconstruction errors, provided by increasing number of iterations, in the case 
of four different binaries (the parameters are indicated in the first column, as explained 
in the text) whose images are generated using the PSF with SR=0.67. As usual MARE 
is a measure of the errors on the magnitudes of the two stars while RMSE is a measure 
of the error on the reconstructed PSF. 





300 it 


4000 it 


8000 it 


12-12 
19 pixels 


MARE 


1.44 x 10~ 3 


1.38 x 10~ 4 


1.31 x 10" 4 


RMSE 


1.8 % 


0.17 % 


0.15 % 


12-16 
19 pixels 


MARE 


5.10 x 10" 4 


1.01 x 10~ 3 


1.52 x 10~ 3 


RMSE 


0.99 % 


0.20 % 


0.27 % 


12-12 
10 pixels 


MARE 


1.26 x 10~ 3 


1.42 x 10" 4 


1.34 x 10" 4 


RMSE 


1.7% 


0.18 % 


0.17 % 


12-16 
10 pixels 


MARE 


1.31 x 10~ 2 


1.92 x 10~ 2 


2.22 x 10~ 2 


RMSE 


1.1 % 


1.3 % 


1.5 % 



values, then our minima were certainly local. 

Before considering other examples it is important to remark that, in the case of the 
binary and of the small star cluster, the stars are reconstructed as single pixels with a 
sufficiently accurate flux value, but the reconstructed images contain artifacts, in the sense 
that other pixels take non zero values. These values are small but they can be disturbing 
in the case of an accurate photometric analysis, for instance of a star cluster, because 
they could be detected as faint stars. Indeed the difference of magnitude between the 
brightest artifact and the true stars is of the order of Am = 8. For this reason we perform 
an analysis of this problem in the case of the binary, using the PSF with the highest SR, 
namely SR=0.67. 

4-3. Artifacts analysis 

We first consider the inverse crime reconstruction of the binary with magnitudes 12-12 
and angular distance 19 pixels. We deconvolve its image generated by means of the PSF 
with SR=0.67 using the same PSF (inverse crime). The algorithm is the standard SGP 
with non-negativity constraint. Also in this case the reconstruction is not free of artifacts 
but they are randomly distributed and their values are quite small: the brightest artifact 
has a magnitude m = 24, hence with a difference Am = 12 with respect to the stars of 
the binary. 

Next we apply the blind algorithm, using as a constraint the exact value of SR and 
we analyze the results provided by the 8000 outer iterations already considered in the 
previous section. In the first two rows of Fig. 3 we show the reconstructed PSF and the 
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MARE = 8.78 x 10~ 5 MAi?E = 4.35 x 10~ 4 MARE = 4.49 x 10~ 4 
# = 120 # = 5 # = 



Figure 3. Binary with magnitudes 12-12 and angular separation of 19 pixels. First and 
second row: the reconstructed PSF and object after 300 (left), 4000 (middle), and 8000 
(right) outer iterations with the exact value of SR (SR=0.67) as a constraint. Third and 
fourth rows: the reconstructed PSF and object after 300 (left), 500 (middle), and 2000 
(right) outer iterations with the underestimated value of SR (SR=0.64) as a constraint. 
The symbol # denotes the number of artifacts. 
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reconstructed object after 300, 4000 and 8000 outer iterations. The object is represented 
using a logarithmic and gray-level reversed scale for stressing the artifacts. After 300 
iterations a few artifacts appear in the reconstructed PSF in a region corresponding to 
the positions of the two stars. Therefore, even if the reconstruction error is small, it 
is evident that it may be convenient to use a larger number of iterations. After 4000 
iterations the PSF artifacts disappear and the reconstruction error is really very small, of 
the order of 0.2%. The situation does not significantly change if we further increase the 
number of iterations. 

As concerns the reconstructed binary, after 300 iterations the artifacts are 
concentrated along two arcs positioned around the two stars. We checked that this 
behaviour is stable if we change the noise realization. Again, if we increase the number 
of iterations the results are better, in the sense that the artifacts are more randomly 
distributed and the intensity of the brightest one decreases. Indeed, after 4000 iterations 
we have Am = 10 and after 8000 iterations Am = 11. We find similar results in the 
case of the open star cluster and therefore we can conclude that a very large number of 
iterations may be required for obtaining very good results, at least if the exact value of 
SR is known. 

However, as briefly discussed in the Introduction, it is not possible to know exactly 
the value of SR. According to astronomers the expected error is of the order of 4 %. 
Therefore, we consider a variation of the constraint of this order of magnitude for images 
generated in the case of SR=0.67; more precisely we consider two values, SR=0.7 and 
SR=0.64. We apply the blind algorithm using as a constraint the corresponding values 
of s. In the first case the reconstructions are definitely worse, the number of artifacts 
considerably increases as well as the error on the PSF. For instance, in the case 12-12 the 
RMSE is of the order of 5 % and does not decrease with increasing number of iterations 
(remember that, in the case of exact value, the error after 8000 iterations is of the order of 
0.15 %). On the other hand, if we underestimate the SR, i.e. we take as a constraint the 
value of s corresponding to SR=0.64, then the results are satisfactory. The reconstruction 
errors for the four binaries already considered in the previous subsection are reported in 
Table [31 By comparing with the results reported in Table [2] and referring to the exact 
constraint, we can conclude that the reconstruction errors are not significantly greater 
than those obtained in the exact case and that the convergence is faster. 

In the case of the binary with magnitudes 12-12 and angular distance 19 pixels we 
show the reconstructions of the PSF and of the binary after 300, 500 and 2000 iterations 
respectively in the third and fourth row of Fig. 3. Quite surprisingly, the artifacts in the 
reconstruction of the binary completely disappear after 2000 iterations and the error on 
the reconstructed PSF is of the order of 2 %. We can add that the same result is obtained 
in the case of the open star cluster. 

In order to further investigate the effect of a wrong value of s, in the case of the 
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Table 3. Reconstruction errors in the case of the four different binaries of Table 2 
(described in the text), considering an underestimated constraint of the blind algorithm 
(SR=0.64). As usual MARE is a measure of the errors on the magnitudes of the two 
stars while RMSE is a measure of the error on the reconstructed PSF. 





300 it 


500 it 


2000 it 


12-12 
19 pixels 


MARE 


8.78 x 10~ 5 


4.35 x 10" 4 


4.49 x 10" 4 


RMSE 


2.1 % 


2.0 % 


2.0% 


12-16 
19 pixels 


MARE 


4.72 x 10~ 4 


3.21 x 10~ 3 


1.62 x 10" 1 


RMSE 


2.0 % 


2.1 % 


4.0% 


12-12 
10 pixels 


MARE 


4.36 x 10" 4 


2.28 x 10" 4 


4.36 x 10" 4 


RMSE 


2.4 % 


2.0 % 


2.0 % 


12-16 
10 pixels 


MARE 


1.99 x 10~ 2 


5.91 x 10~ 2 




RMSE 


2.4 % 


3.5 % 


4.0 % 



binary 12-12 and SR=0.67, we also consider the cases SR=0.4, 0.6, 0.8 and 1. The error 
on the PSF is of the order of 22 % in the case SR=0.4, about 5 % in the case SR=0.6 and 
about 11 % in the two other cases. We can add that the artifacts disappear in the case of 
underestimated SR. In conclusion an underestimate of SR of about 10 % (which does not 
correspond to the precision achievable in the experimental estimation of this parameter) 
is still acceptable, while an overestimate can be dangerous in all cases. 

4-4- Complex and diffuse objects 

As additional examples of astronomical targets, we consider three HST images: the Crab 
nebula NGC 1952, the galaxy NGC 6946 and the planetary nebula NGC 7027. In all 
cases we assume an integrated magnitude equal to 10 and, for each one, we obtain three 
blurred images by convolving with the three PSFs of SR=0.67, 0.40 and 0.17. Again a 
background in H-band is added to all images and the results are perturbed with Poisson 
and additive Gaussian noise. In the first column of Fig. 4 we show the three objects in 
reversed scale of gray levels while in the second column we show their blurred images in 
the case SR=0.67. 

As concerns the initialization of the blind algorithm we use the same PSFs already 
used in the case of stellar objects, namely obtained by suitable autocorrelations of the 
ideal PSF of the telescope. However in the case of complex and diffuse objects, as already 
remarked by other authors (see, for instance, pH]), a difficult and crucial point is the 
choice of the number of inner iterations. We do not have a rule which can be successfully 
applied to all cases as for the stellar objects, i.e. (nf,n h ) = (50, 1). By several attempts 
we find "best" numbers for each case, but it is obvious that this is not a satisfactory 
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Table 4. Reconstruction errors for complex and diffuse objects. In third and fourth 
columns, the best errors achieved by SGP with the true and the initial PSFs, respectively. 
In the fifth column, the best error obtained using a maximum of 100 outer iterations (for 
the choice of the inner iterations see the text). Finally, in the last two columns, the error 
between the true PSF and the initial one, followed by the error between the true PSF 
and the one obtained in conjunction with the best reconstruction of the corresponding 
object. 



Image 


SR 


RMSE obj 


RMSE° bJ 


RMSE° bJ 


RMSE[ sf 


RMSE% S} 




0.67 


11 % 


16 % 


12 % 


32 % 


6.7% 


Crab 


0.40 


12 % 


20 % 


14 % 


54 % 


12 % 




0.17 


15 % 


22 % 


16 % 


55 % 


16 % 




0.67 


14 % 


23 % 


16 % 


32 % 


7.4 % 


Galaxy 


0.40 


16 % 


30 % 


19 % 


54 % 


16 % 




0.17 


20 % 


35 % 


23 % 


55 % 


21 % 




0.67 


3.2 % 


6.8 % 


6.8 % 


32 % 


32 % 


Nebula 


0.40 


3.5 % 


8.9 % 


8.9 % 


54 % 


53 % 




0.17 


4.2 % 


9.0 % 


7.9 % 


55 % 


43 % 



situation. For instance, in the case of the Crab nebula we find (nf,nh) = (13,22) for 
SR=0.67 and (n/,n/i) = (13,27) for SR=0.40, if we search for minimum RMSE on the 
object using 100 outer iterations. Moreover, even if the algorithm is convergent, the limit 
is not, in general, a sensible solution: a suitable stopping of the outer iterations is required. 
In other words the algorithm is semi-convergent [341 H] as concerns the outer iterations, 
i.e. the RMSE on the object first decreases, reaches a minimum and then goes away. We 
do not have a proof of this feature which derives from the numerical experiments. It is 
obvious that such a situation is not satisfactory and we briefly discuss this point in the 
next section (see also the Introduction). 

In Table H] we report the best results we have obtained while in the third column of 
Fig. 4 we show the reconstructions of the objects provided by the blind algorithm. We 
stress the case of the planetary nebula: it seems that in this case the algorithm is unable 
to improve the reconstruction with respect to that provided by the initial guess. We also 
remark that the error on the reconstructed PSF depends on the object: for instance, for 
the galaxy it is larger than that for the Crab. 

We conclude by reporting an experiment intended to check the effect of an 
underestimated or overestimated SR in the reconstruction of a diffuse object. We consider 
the case of the Crab nebula and PSF with SR=0.67, and we apply our algorithm with 
SR=0.6 and SR=0.8. In the first case, the minimum RMSE on the object and the PSF are 
equal to 13% and 8.2%, while in the other the errors on object and PSF are 12% and 6.7%, 
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Figure 4. First column: the objects used for image generation; second column: the 
blurred images in the case of SR=0.67; third column the reconstructed objects obtained 
with our blind algorithm. First row: the Crab nebula NGC 1952, second row: the spiral 
galaxy NGC 6946, third row: the planetary nebula NGC 7027. 



which are essentially the same values obtained with the correct SR. The small variance 
observed suggests that, in presence of complex objects, the availability of a correct SR 
does not represent a crucial point. In this case, an overestimate of the SR value seems to 
be preferable. 

5. Concluding remarks and perspectives 

In this paper we propose a blind algorithm, based on SGP, for the reconstruction of 
astronomical images acquired by a telescope equipped with an AO system. The algorithm 
can be classified as an inexact alternating minimization of the KL divergence depending on 
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both the object and the PSF. The crucial point is the introduction of different constraints 
on the object and PSF and this can be done in a correct way by using the particular 
structure of SGP. Moreover the convergence of the algorithm is derived from general 
results proved in [2] on the convergence of inexact AM. 

Since the problem is non-convex, the limit points of the sequence of iterates may 
depend on several parameters and, more specifically, on the initialization of the outer 
iterations and on the numbers of inner iterations. In the case of stellar (point-wise) 
objects we have rules for the initialization and the numbers of inner iterations which seem 
to be suitable for all cases we have considered. Obviously, the effectiveness of the approach 
must be tested by application to a much broader set of simulated images in a wide program 
of simulations (see, for instance, Chap. 12 of [I]) as well as to real images. We also observe 
that our blind approach can be possibly used in conjunction with specific codes, such as 
the so-called StarFinder [22] , developed for accurate photometric and astrometric analysis 
of star clusters. These codes contain a method for extracting the PSF from the image of 
the star field; this PSF can be compared with and/or replaced by that provided by our 
blind approach; also in this case the analysis of simulated star fields, in particular crowdy 
star fields, can help to understand when the blind approach is required. 

As already remarked in the Introduction, the situation is not so clear in the case 
of more complex astronomical targets. Our preliminary simulations indicate that the 
proposed blind method has the semi-convergent property as concerns the outer iterations 
(the numbers of the inner iterations are fixed by the user and, in any case, they should 
not be too large). The difficulties found in optimizing the numbers of inner iterations 
may reside in the fact that in this paper we do not introduce regularization in the 
objective function. Due to the flexibility of SGP the method can be easily generalized to 
differentiable regularizers (one only needs to compute the gradient of the penalty and its 
positive part) and this generalization will be the subject of future work. We stress again 
that the crucial point is to identify regularizers which are suitable for specific classes of 
astronomical targets as well as regularizers which are suitable for AO corrected PSFs. 

The codes of the algorithms presented and used in this paper are available under 
request. 
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